/*
 * file: benchmark.c
 * author: adapted from CSE6230
 * desc: benchmark performance of openmp naive matrix multiplication
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <float.h>
#include <math.h>
#include "matrix_utils.h"
#include "matmult.h"
#include "stopwatch.h"
#include "mkl.h"

#define GIGA 1.0e9
#define N_TEST (sizeof(test_sizes)/sizeof(int))
#define N_RUNS 1

/* init some vars *****************************************************/
const char trans = 'N';
const double alpha = 1.0;
const double beta = 1.0;
static struct stopwatch_t *timer;
static double *A1, *A2, *B1, *B2, *C1, *C2;

/* forward decls *****************************************************/
void benchmark_init(const int);
bool benchmark_val(const int, const int);
double benchmark_time(const int);
void benchmark_clean();

/* main **************************************************************/
int main(int argc, char *argv[])
{
    //1) init stuff
    int test, pass, i;
    int test_sizes[argc-1];
    for (i=1; i<argc; i++) { //load comm args
        test_sizes[i-1] = atoi(argv[i]);
    }
    int max = test_sizes[argc-2];
    double perf_total = 0.0;
    double avg;
    benchmark_init(max);

    //2) validate implementation
    printf("\n**Validating**\n");
    
    for (test=0; test<(int)N_TEST; test++) {
        int n = test_sizes[test];
        if (benchmark_val(max,n) == false) {
            printf("\tValidating n=%d.....fail\n",n);
        } else {
            printf("\tValidating n=%d.....pass\n",n);
        }
    }

    //3) time implementation
    printf("\n**Timing**\n");
    for (pass=0; pass<(int)N_RUNS; pass++) {
        for (test=0; test<(int)N_TEST; test++) {
            int n = test_sizes[test];
            double gflops = benchmark_time(n) / GIGA;
            printf("\tn=%d: %f gflops/s\n",n,gflops);
            perf_total += gflops;
        }
        printf("\n");
    }
    avg = (perf_total/N_TEST)/N_RUNS;
    printf("\n**Average performance: %f gflops/s**\n", avg);

    //4) cleanup
    benchmark_clean();
    return 0;
}

/* benchmarking steps *************************************************/
void benchmark_init(const int max) 
{
    //allocate matrices
    A1 = (double *)malloc(sizeof(double)*max*max);
    A2 = (double *)malloc(sizeof(double)*max*max);
    B1 = (double *)malloc(sizeof(double)*max*max);
    B2 = (double *)malloc(sizeof(double)*max*max);
    C1 = (double *)malloc(sizeof(double)*max*max);
    C2 = (double *)malloc(sizeof(double)*max*max);
    if (A1 == NULL || B1 == NULL || A2 == NULL || B2 == NULL || C1 == NULL || C2 == NULL) {
        printf("Error: Allocation failed for at least one matrix\n");
        exit(-1);
    }
    
    //start timer
    stopwatch_init();
    timer = stopwatch_create();
}

bool benchmark_val(const int max, const int n) 
{
    //generate a random matrix of size n, and verify our implementatino
    //against MKL's dgemm
    random_matrix(n,n,A1);
    random_matrix(n,n,B1);
    copy_matrix(n,n,A1,A2);
    copy_matrix(n,n,B1,B2);
    zero_matrix(n,n,C1);
    zero_matrix(n,n,C2);
    dgemm(&trans,&trans,&n,&n,&n,&alpha,A1,&n,B1,&n,&beta,C1,&n); //mkl
    matmult(n, A2, B2, C2);
    
    return verify_matrix(n,n,C1,C2);
}

double benchmark_time(const int n)
{
    //init
    int i;
    long double flops, t;
    long double ops = (long double)2.0*(long double)n*(long double)n*(long double)n;

    random_matrix(n,n,A1);
    random_matrix(n,n,B1);
    random_matrix(n,n,A2);
    random_matrix(n,n,B2);
    zero_matrix(n,n,C1);
    zero_matrix(n,n,C2);

    //burn in with 2 runs
    for (i=0; i<2; i++) {
        matmult(n, A1, B1, C1);
        matmult(n, A2, B2, C2);
    }

    //start and stop timer around matmult calls
    stopwatch_start(timer);
    for (i=0; i<N_RUNS; i++) {
        //alternating matrices to clear the cache
        matmult(n, A1, B1, C1);
        matmult(n, A2, B2, C2);
    }
    stopwatch_stop(timer);

    //return results
    t = stopwatch_elapsed(timer) / (long double) N_RUNS;
    flops = (long double) ops / t;
    return (double) flops;
}

void benchmark_clean()
{
    //free resources
    free(A1);
    free(A2);
    free(B1);
    free(B2);
    free(C1);
    free(C2);
    stopwatch_destroy(timer);
}
